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Abstract 



We study the phase diagram of the ferromagnetic g-state Potts model 
on the various three-dimensional lattices for integer and non-integer values 
of g > 1. Our approach is based on a thermodynamically self-consistent 
Ornstein-Zernike approximation for the two-point correlation functions. We 
calculate the transition temperatures and, when the transition is first order, 
the jump discontinuities in the magnetization and the internal energy, as well 
as the coordinates of the critical endpoint in external field. Our predictions 
are in very good agreement with best available estimates. From the numerical 
study of the region of weak first-order transition, we estimate the critical value 
Qc for which the transition changes from second to first-order. The g — > 1"*" 
limit that describes the bond-percolation problem is also investigated. 
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I. INTRODUCTION 



The Potts model is a generalization of the Ising model to more (or less) than two 
states per spin. This model is realized in many different experimental situations (see Ref. 
for a review) and it has been the subject of considerable theoretical attention over the 
last two decades. At low enough temperature and in the absence of a symmetry breaking 
magnetic field, the ferromagnetic Potts model undergoes a transition from a disordered 
phase in which all states are equally populated to an ordered phase that is characterized 
by a non-zero spontaneous magnetization associated with preferential occupation of one of 
the q states. The order of the transition and the critical properties crucially depend on 
the spatial dimension d and the number of states q. Whereas the mean-field description 
predicts the transition to be first-order for all g > 2 irrespective of the dimension 0, the 
exact analysis in two dimensions shows that the transition is continuous for g < 4 [Q. No 
exact results are available in three dimensions but there is by now general consensus that the 
3-state Potts model on the simple cubic lattice undergoes a weak first-order transition with 
a very small latent heat and a large correlation length 0. Since the transition is continuous 
for g < 2 (g = 2 is the standard Ising model whereas the limit g — > 1"*" corresponds to 
the bond-percolation problem), the order of the transition in 3-d must change at some gc 
between 2 and 3. There are, however, only few estimates of gc, and the behavior in the 
neighborhood of gc is not precisely known 

To gain a better perspective on this issue and, more generally, to determine accurately 
the phase diagram of the ferromagnetic Potts model on the various three-dimensional lat- 
tices, we propose an analytical treatment that allows us to investigate integer as well as 
non- integer values of g (g > 1). Our approach derives from an approximation that was first 
introduced to study the phase behavior of simple fluids and that has been successfully 
applied to a variety of spin models, including the spin-i Ising model fiy], the Blume-Capel 



model ||TT| , and several models of random magnets |T^. This approximation, called the 
" self-consistent Ornstein-Zernike approximation" (SCOZA), is based on what is known in 
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liquid-state statistical mechanics as the Ornstein-Zernike formalism, namely (i) the exact 
Ornstein-Zernike equations that relate the two-point connected correlation functions to the 
two-point direct correlation functions (the former are obtained as functional derivatives of 
the field- dependent free-energy and the latter as functional derivatives of the Legendre- 
transformed magnetization-dependent Gibbs potential) and (ii) an Ornstein-Zernike ansatz 
that assumes that the direct correlation functions have the same range as the interaction 
pair potential. The dependence of the correlation functions on the control variables (tem- 
perature and magnetization) is determined by two types of constraints: on one hand, exact 
zero-separation conditions for the connected correlation functions; on the other hand, ther- 
modynamic self-consistency which ensures that the same Gibbs free energy is obtained from 
the two-point correlation functions, whether one integrates the internal energy with respect 
to the inverse temperature or one integrates twice the inverse susceptibility with respect to 
magnetization. 

The general SCOZA formalism for the ferromagnetic g-state Potts model in the presence 
of g — 1 external magnetic fields is presented in section In section III, we restrict attention 
to the case of a single symmetry breaking field that favors one state while preserving the 
symmetry between the q — 1 remaining states. This allows us to study the thermodynamics 
of the model for any real value of g > 1 by solving a partial differential equation with 
temperature and magnetization as independent variables. The numerical results for the 



various 3-d lattices are presented and discussed in section We first study the cases g = 3 
and g = 4 and compare to the best available estimates. We next consider 2 < g < 3 and 
provide a numerical estimate of the critical value Qc at which the transition changes from 
second to first order. Finally, the case 1 < g < 2 is discussed and the limit g 1+ is 
investigated. 

II. SELF-CONSISTENT ORNSTEIN-ZERNIKE APPROXIMATION FOR THE 

FERROMAGNETIC POTTS MODEL 
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A. Definitions and exact relations 



The gr-state Potts model on a d-dimensional lattice of N sites is defined by the general 
Hamiltonian 

N q-1 

n^-jJ2 i<i^s.,s, - 1) - E E ^"(«^^»." - 1) (1) 

<ij> i=l a=l 

where J is a positive coupling parameter (we consider here only ferromagnetic interactions) , 
5siSj is the Kronecker symbol, and the sum on < ij > runs over distinct nearest-neighbor 
(n.n.) pairs. Each site variable Sj can assume q possible values a — 1 . . .q and the /if s are 
site-dependent magnetic fields that couple to the distinct states. (Because of the condition 
Sa=i ^si,a = 1, only q — 1 fields are relevant and one may take hj = without loss of gener- 
ality.) These site-dependent fields have been introduced in order to generate the correlation 
functions (see below), but ultimately we are interested in the system in the presence of a 
uniform field that favors only one state. The field-dependent free enegy JF is then defined 

by 

= jr(/? J, {h^}) = - In Tr exp[-/??i] (2) 

where /3 = 1/(A;bT) is the inverse temperature and = (3hf. The derivatives of J-' with 
respect to the external fields generate the local magnetizations mf, 

dT 

mr=<ar>=--^, a = l,2,...,(g-l) (3) 

and the connected spin-spin correlation functions 



Gtf =< <<^^ > - < < X -^f-^ a,P=l,2,...,{q-l) (4) 

where wc have introduced the spin variables = qSg^^a — 1- Because of the condition 
Yla=i'^i — 0, there are only (g — 1) independent magnetizations (per site) and q{q — l)/2 
independent correlation functions {C^f — G^f ) with 

q-1 



a=l 



and 

G^' = -E< a = l,2,...,^. (6) 

0=1 

Since we are interested in working with the magnetizations rather than the fields as 
control parameters, we introduce the Legendre transformed Gibbs potential 

9-1 
i a=l 

which generates the direct correlation functions C^"^, 

Cff^ a,/3= l,2,...,(g-l). (8) 

The two-point correlation functions and the direct correlation functions C^^ are related 
by a set of Ornstein-Zernike (OZ) equations that result from the Legendre transform, 

q-1 

E E ^"^k^ = '^".A.- /3 = 1, 2, {q-l). (9) 

A; 7=1 

In the limit of uniform fields, one has = and the two-point correlation functions 
depend only on the vector r that connects the two sites. The OZ equations then take a 
simple form in Fourier space, 

q-1 

J2 C"^(k)G^''(k) = Sa,p 2, {q-1), (10) 

7=1 

and Eq. (8) yields the (inverse) susceptibility sum-rules: 

^°'(''='') = S «.,3=1.2....,(,-1). (11) 

Owing to the properties of the spin variables, the values of the correlation functions 
G°-^{v) at r = can be expressed in terms of the magnetizations as 

G'"''(r = 0) = -(l + m")(l + m^) a ^ (3 (12) 

^""(r^O) = -(l + m")(l + m"-g). (13) 
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On the other hand, the enthalpy of the system (which is just the internal energy in the ab- 
sence of external field) depends on the values of the correlation functions at n.n. separation. 
Using the identity 6s,s, = ZlLi ^s.a^s.a, we get 

^ = H = E < (^'^-.«. - 1) >= = e) + {m-Y] , (14) 

^ <ij> ^ a=l 

where c is the coordination number of the lattice, and e denotes a vector from the origin to 
one of its nearest neighbors. 

For the theory to be thermodynamically consistent, the Gibbs potential Q must have the 
same value when obtained via the susceptibility route, i.e., Eq. (11), or via the energy (or 



enthalpy) route, i.e., Eq. (|1J)- This thermodynamic consistency is embodied in Maxwell re- 
lations between the partial derivatives of Q with respect to the independent control variables. 
In the present case, these relations have the following form: 

^C-^(k = 0) = -l^^tlG-'-'(r = e) + (m^^| a. /J = 1. 2. ....(,- 1) . (15) 

7=1 

which derive from Eqs. (11) and (14) (from now on the coordination number c is absorbed 
in the inverse temperature variable A = (3cJ). These partial differential equations (PDE) 
extend the self-consistent equation considered in Ref. 0] for the spin-| Ising model to (g — 1) 
order parameters. 

The theory presented in this paper is based on approximations for the structure of the 
direct correlation functions and the above exact equations are the starting point of our 
analysis. 



B. Ornstein-Zernike approximation 

In general, the direct correlation functions are expected to remain of finite range, even 
at a critical point |]T3[. Following the Ornstein-Zernike approximation, we assume here that 
they have exactly the range of the exchange interaction in the Potts Hamiltonian, Eq. (1), 
which implies that they are truncated at nearest-neighbor separation. We thus write 
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C-P{r) = cf5,,^ + Cf5,,,, (16) 

or in Fourier space, 

C"'^(k) = Cf + c Cf^A(k) (17) 



where A(k) = 'Y^^e^^'^ is the characteristic function of the lattice. The q{q — 1) distinct 
functions Cq^ and C^^ depend on the inverse temperature and on the magnetizations, but 
this dependence is not given a priori. As is well-known in liquid-state theory |[T^, when 



one assumes some approximate but explicit dependence of the direct correlation functions 
upon the thermodynamic variables, as in the random phase approximation or in the mean- 
spherical approximation, the theory is in general thermodynamically inconsistent. In our 
self-consistent OZ approximation, the functions Cq^ and Cf^ are determined by enforcing 
the zero-separation conditions, Eqs. (0)-(|l3), and by imposing thermodynamic consistency 



via Eq. ([T^). 

The above equations, Eqs. (10)-(17), define the SCOZA for the Potts model in the most 
general situation, i.e., in the presence of g — 1 independent ordering fields. To make it a 
workable scheme, one has to ensure that the number of equations equals the number of 
unknowns and that the PDE's are supplemented with proper boundary conditions. This is 
the case for q = 2 (the Ising model), and more generally for integer values of g > 2. Indeed, 
Eqs. ([12|) , (p!3|), and ([15|) provide q{q — 1) independent equations for the q{q — 1) unknowns, 
and the boundary conditions for the g-state model may be obtained either exactly (for A = 
and m° = g — 1) or from the solution of the (g — l)-state model (for = — 1). Therefore, 
when g is an integer, one can solve the equations recursively. In practice, however, this is a 
difficult numerical problem, as one has to solve g(g — 1)/2 coupled PDE's with g independent 
variables, the g — 1 independent magnetizations plus the inverse temperature (the equations 
must be solved in the full g-dimensional space even if one is interested ultimately in the zero- 
field case only). It is thus necessary to develop a simpler approach in which the number of 
PDE's does not depend on g. Moreover, since our goal is to study the Potts model for integer 
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as well as non-integer values of g, one should be able to perform an analytic continuation of 
the number of states to arbitrary real values. 



III. A SCOZA FOR NON-INTEGER VALUES OF q 

A. The partial differential equation 

To be able to study non-integer as well as integer values of g, we consider the case of 
a magnetic field that acts only on one state, say state 1, while preserving the symmetry 
between the remaining (g — 1) states. Without loss of generality, we take = h and 
h? = = . . .h'^ = Q. It follows that m? = = . . . = m'^ = —m}/ (g — 1), and there is only 
one independent order parameter, m} = m. By symmetry, one has also only four distinct 
connected (resp. direct ) correlation functions, for instance G^^, G^^, G^^, and G^^ (resp. 
(jii^ (j22^ (ji2^ g^^^ C^'^). It is then easy to block-diagonalize the matrices of correlation 
functions so to rewrite the OZ equations, Eq. ([l^), as 

/ 



^ G"(k) G^\k) ^ 

^{q-2)G'\k) [G22(k) + (g-3)G^3(k)]^ 



C"(k) Ci2(k) 
(g-2)ai2(k) [d''ik) + iq-3)C''ik)] 



(18) 



[G^\k) - G^%k)] [d^\k) - C^^ik)] = 1. (19) 
Moreover, Eq. (6) implies that 

^11 
g — 1 

^''=(g-l)(g-2)"^' ^^^^ 
from which we obtain, by replacing in Eqs. ( [T^ ) and (|T^, 

G^^ = G^^ = ]^G^\ (22) 
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We finally get 

G'\k) = ^—^ , (23) 

and 

e-(k)= «i^ + 2i^^. (24) 

The enthalpy S (see Eq. (|T^)) is now given by 

S = -^[G'\r = e) + (g - l)G'\r = e) + -^m% (25) 
and since there is only one order parameter, there is only one inverse-susceptibility relation, 

The partial differential equation expressing thermodynamic consistency follows from 
Eq. (0): 

|P"(k^O)-li^e-(k^O)]^-^ (27) 

^ ^'.[G'"(r = e) + (g-l)G22(r = e)]. 



2q dm? 

So far, all equations are exact. 

The OZ approximation, Eq. (|1^, introduces 4 unknown functions, Cq^, Cq^, C\^, Cp, 

C"(k) = Co^i + cC|^A(k) (28) 
522(k) =C22 + cCfA(k) . (29) 



Introducing the variables 



"'^^11 _ 122:2^22 V-^^J 
'-'0 2 5-1^0 



and 



22 



(72 

= , (31) 
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we have 



2q-l 



2q-l 



and 



(32) 



(33) 



Using Eqs. p3|j2^) and introducing the lattice Green's function |jT4 

1 -ik.r 

Piz^r) = — ^ / d^'k^^ , 



(2vr) 



(34) 



1 - z\{k) 

we finally obtain the expressions of two-point connected correlation functions in real space, 

P(z.r) 



r^U _ 1 9-2^22 
^0 2g-1^0 



(35) 



and 



G^^{r) q~2P{z\r) 



+ 2 



(36) 



(g_l)2 ■ C'22 • 

The two unknown quantities, Cq^ and Cq^, are eliminated by enforcing the exact zero 
separation conditions following from Eqs. (12) and (13), namely 



^^^(r = 0) = -(1 + m)(l + m - g). 



(37) 



This yields 



[{q — 1)^ + m] [1 + m — q] 



G^'{r) = (1 +m)(g- 1 



m] 



P(z) 



(38) 



(39) 



G22(r) 



G^\r) ^q{q-2){q-l-m)P{z*,r) 



(40) 



where P{z) 
obtained: 



{q-iy (g-1)' Piz") ' 

P{z, r = 0). By inserting the above results in Eq. (0), the following PDE is 
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A|P(,)(i_,))._(l±i!|l_i^ (41) 

A&^^^o-'- '""'^ + w ^ - ^' wfi" 

where the relation P(z, r = e) = {P{z) — l)/z has been used. We see that q is now a 
parameter that can take any real value. Setting g = 2 in Eq. ( PT| ) gives back the PDE for 
the Ising model with only one unknown function, z{\^ m) |jlO|]. In the general case, however, 
there are two unknown functions, 2;(A, m) and z*{\, m), but only one equation relating them. 
Additional input is therefore needed. We address this point as well as the required boundary 
conditions for Eq. (^) in the following. 

B. Boundary conditions and further approximation 

The full range of variation of A is from to +oo and that of the magnetization m 
goes from —1 to (g — 1). The initial condition for Eq. pT]) at A = is provided by the 
exact solution of the noninteracting Potts model in an external field. One has in particular 
Cl\X = 0, m) = CP(A = 0, m) = 0, so that z{X = 0, m) = z*{X = 0, m) = 0. Similarly, the 
condition for m = q — 1 follows from the fact that all site variables are in state 1 {h oo), 
which again leads to z{X, m = q — 1) = -2*(A, m = q—1) = because Cq^ and Cq^ +oo. At 
the boundary m = — 1, the state 1 is suppressed: the system then corresponds to a (g — 1)- 
state Potts model in zero field, whose solution is unknown in general. However, in order to 
study the zero-field transition of the g-state Potts model (g > 1) that leads to preferential 
occupation of state 1, it is sufficient to consider the region m > 0. The boundary condition 
at m = results from the Zq symmetry of the Hamiltonian in Eq. when the external 
fields are turned off: above the transition temperature, one must have h{X,m = 0) = and 
all g states must be equivalent. As a consequence, the correlation functions G^^ and G^^, 
and similarly C^^ and C^^, are equal; hence, Cq^(A, m = 0) = Co^(A, m = 0) and Cl^{X, m = 
0) = CP(A,m = 0), which implies via Eqs. (^) and (§T]) that z{X,m = 0) = z*{X,m = 0). 
The condition h{X^m = 0) = can be reexpressed as a constraint on the enthalpy £. 
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Indeed, since the boundary condition at A = garantees that h{\ = 0, m = 0) = 0, and 
since thermodynamic consistency imphes that dh/d\\m = l/{NcJ) dS/dm\x, the zero-field 
condition in the disordered phase requires that dS /d'm\x,m=o = 0. 

Having settled the boundary-condition problem, we must still find an additional equation 
that, together with Eq. (0), allows to uniquely determine the functions z{X, m) and z*{X, m). 
We do so by considering the exact expressions of the direct correlation functions at nearest 
neighbor separation, and Cp, at the two boundaries m = and m = q — 1. We already 
saw that C|^(A,m = 0) = Cf^{X,m = 0) as a result of the Zq symmetry. The calculation 
for m = g — 1 is detailed in Appendix A and leads to the following expression: 



ACiX) = Cl\X,m = g - 1) - Cf (A,m = q - I) = -^{e^"'" - !)'• (42) 

We then make a linear interpolation between the m = and m = q — \ values for all 
temperatures: 

Tf) 

Cl\X,m) - Cl\X,m) ^ ^--^AC(A) , (43) 
which yields the following additional equation for z(A,m) and z*(A,m): 

777 

zPiz) - (1 + m^z^Piz") = -c(l + m)(g - 1 - m)- rAC(A), (44) 

(9-1) 



where we have used Eqs. (|35D-(|38D to eliminate C^^ and Cq^. 

C. Numerical procedure 

The numerical integration of Eq. (|4l|) , coupled with Eq. (0), is carried out by an explicit 
algorithm, the partial derivatives being approximated by finite difference representations. 
Eq. is rewritten as a nonlinear diffusionlike equation, 

dz (l + m)(g— 1 — m) Id'^AS , . 

9A " ~ (g - 1)[(1 - z)P'{z) - P{z)] ^ ^2~d^' ^ ' 

with P\z) = dP{z)/dz and AS = {q - 1 - m){{l + m)[P{z) - l]/[zP{z)] + {q-2)[P{z*) - 
l]/[z*P{z*)]}. The derivative of z with respect to A is obtained from 
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where SX is the elementary step in the "time-hke" direction, and the second derivative of 
AS with respect to m is given by 



d^AS 



dm? 



A£{X, m + Sm) - 2A£{X, m) + A£{X, m - 5m) 

1 y^') 



where 6m is the elementary step in the "space-like" direction. 

One starts from A = for which z = z* = AS = (noninteracting system). Then, 
Eq. (|45| ) gives z at the next step in the A-direction, A = 6X. At this temperature, z* is calcu- 
lated from Eq. (0) by a Newton-Raphson algorithm. Once z and z* are known, AS and its 
second derivative are computed. At the next step, the temperature is decreased (A X + 6X) 
and the same procedure is repeated. 6X is gradually decreased as the spinodal is approached 
(see below) to ensure numerical stability of the explicit scheme. As previously mentioned, the 
boundary condition at m = g — 1 is known exactly, with z{m = q — 1) = z*{m = q — 1) = 0, 
and the zero-field condition at m = is enforced by setting AS\m=Sm = AS\m=-6m which 
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readily yields from Eq. (47) d'^AS/dm^\^^^ = 2[AS{X,6m) - AS{X,0)]/6m 

It follows from Eqs.(26) and (32) that the inverse susceptibility goes to zero when z —>■ 1. 
(Note that from the definition of the lattice Green's function, Eq. (p^), z and z* must stay 
in the interval [0, 1].) The locus of the points for which z = 1 defines a spinodal curve in the 
m — T plane. We denote Tc the temperature at which the spinodal has its maximum. When 
the maximum of the spinodal occurs at m = 0, the zero-field transition is continuous and the 
spontaneous magnetization, obtained from the condition h{m,T) = d{Q / N) / dm = 0, varies 
with T as depicted schematically in Fig. (la). When the maximum of the spinodal occurs 
at m = m^ 7^ 0, the zero-field transition is first-order and the spontaneous magnetization 
jumps to a non-zero value Am at a transition temperature Tt < Tc, as shown in Fig. (lb). 
The transition remains first-order in the presence of a small non-zero external field and for 
h = he it ends in a second order critical point. This critical endpoint corresponds to the 
maximum of the spinodal. 
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In the numerical procedure, the spinodal curve is never reached exactly. With the choice 
of a small parameter e (typically e = 10~^ or 10~^), we define a pseudo-spinodal as the 
points of the grid formed by the discretized values of m and A that satisfy 1 < z < 1 — e. 
The pseudo-spinodal then serves as a boundary condition for the finite-difference equation 
in the low-temperature phase. In this phase, 6X is fixed to a very small value (typically 
10^® or 10^^). We have checked that the location of the spinodal does not move when e is 
changed from 10~^ to 10~^ and that the whole calculation is stable when 6X is decreased 
from 10~^ to 10^''. All calculations are done with 6m = {q — l)/200. 

The numerical integration has been carried out for the sc, bcc, and fee lattices. The 



corresponding Green's functions P{z) are well documented |jT^ and can be easily tabulated. 
To test the numerical procedure, we have checked that the results for the Ising model 
are recovered with good accuracy. For instance, we have found l3cJ = 0.22126 for the 
inverse critical temperature on the sc lattice, in perfect agreement with the value obtained 
in Ref. 0. 



IV. RESULTS AND DISCUSSION 
A. q = 3 and 4 

We have first considered the 3 and 4-state Potts models that have been much studied in 
the recent literature (the 3-state model plays a important role in both condensed matter and 
high-energy physics, and the 4-state model is a special case of the Ashkin- Teller model). On 
the three lattices, we have clearly found that the system undergoes a first-order transition 
with a critical point in non-zero magnetization (or field). For illustration, we show in 
Fig. 2 the spinodal and spontaneous magnetization curves for the sc lattice. The (inverse) 
transition temperatures, given in Table I, are in very good agreement with the most recent 
estimates from Monte carlo simulations, q(3tJ = 0.550565(10) for the 3-state model on the 
sc lattice §] and q/3tJ = 0.628616(16) and 0.455008(10) for the 4-state model on the sc and 
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bcc lattices, respectively [0. (To the best of our knowledge, the transition temperatures for 
the fee lattice are given here for the first time.) Remarkably, it appears that the precision of 
the theory (0.05%) is better than for the Ising model ||T0[, despite the use of the additional 
approximation for — Cp which has led to Eq. (44). This is probably due to the fact the 
correlation length, although large, remains finite at the transition. 

We also give in Table I the predictions for the dimensionless latent heat Ae = e(Tj^) — 
e{Tf), where e = (c/2) < ds^sj >, and for the jump discontinuity in the order parameter 
at the transition Am/{q — 1) (with our definition of the local magnetizations in Eq. (3), 
m varies here between and q — 1). For the 3-state model on the sc lattice, the predicted 
value of the latent heat is somewhat larger than the estimate from Monte Carlo simulations, 
Ae = 0.1606 ± 0.0006 |T^, but it is smaller than the estimate from series expansions, 
Ae = 0.264 ± 0.011 [1^. (As noted in Ref. [0, the large difference between these two 



estimates is surprising and has received no explanation.) Similarly, our value of Am/2 is in 



the range defined by the Monte Carlo estimate, Am/2 = 0.395 ± 0.005 ||T8|, and the series 



expansion estimate, Am/2 = 0.505 ||T^. Our results also show that Am is only weakly 
dependent on the lattice structure. This quasi-universality is a consequence of the large 
value of the correlation length at the transition. (In two dimensions. Am is the same for 
the square, triangular and honeycomb lattices, but this is a consequence of the star-triangle 
relation pf.) 

For completeness, we also report in Table I the coordinates of the critical endpoint in 
non-zero external field. For the 3-state model on the sc lattice, the theoretical predictions 
are very close to recent Monte Carlo estimates, qPcJ = 0.54938(2) and Pchc = 0.000775(10) 
20| . We may note, however, that the accuracy for the critical temperature is "only" 0. 



0. 



B. 2 < g < 3 

As q is decreased below 3, Am and he decrease and Pc gets close to Pt, as shown in 
Table II in the case of the sc lattice. This indicates that the first-order character of the 
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zero-field transition becomes weaker and weaker. For q < 2.4, the spinodal is extremely 
fiat and to within the accuracy of our calculations we cannot distinguish /3t from /3c nor he 
(or Am) from zero. Therefore, the direct analysis of the numerical results cannot yield a 
clear-cut conclusion regarding the order of the transition in this range of q. On the other 
hand, the asymptotic analysis of the equations in the vicinity of a putative critical point 
indicates that because of the linear interpolation formula for the direct correlation functions 
at n.n. separation, Eq. (43), the transition is first-order for all 5 > 2: indeed, as shown 
in Appendix B, this additional approximation is incompatible with the zero-field condition 
dS / dm\x^rn=o = for g > 2 if the transition is continuous. This latter conclusion is likely 
to be an artifact of the too crude linear interpolation formula. A proper resolution of the 
question of the existence of a critical value > 2 would require a fully self-consistent SCOZA 
that does not rely on an approximate, explicit relation between z and z*. In the absence of 
such a theory, interesting insight into the potential change of the transition from first- to 
second-order can still be gained by studying the numerical results in the region of q where 
the first-order character of the transition can be clearly identified with numerical precision, 
i.e., for q > 2.4. As illustrated by the very good accuracy of the theoretical predictions for 
q — 3, we expect this region to be less dependent on the detail of the relation between z and 
z*, and the putative changeover behavior at qc can then be investigated by extrapolation of 
the data to q less than 2.4. 

This procedure has been used to analyze the data shown in Figs. 3 and 4 for the 
latent heat Ae and the second-moment correlation length ^ in the disordered phase at the 
transition point (^ is defined by G^^(k) ~ G^^(0)(1 + ^'^k'^),k — > 0; specifically, one has 
c^"^ — z/{l — z)). For the sc lattice, we find that Ae is well fitted in the range 2.4 < q < 2.8 
by the simple algebraic form Ae ~ a{q — qc)" with q^. 2.15 and a ~ 1.70 (if the result for 
g = 3 is also taken into account, we find instead qc ~ 2.21 and a ^ 1.44). Similarly, the 
rapid increase of ^ is well described by ^ ~ b{q — qc)~^ with qc ~ 2.14 and u 1.14 (this 
fit is not modified when the result for g = 3 is added, as can be seen in Fig. 4). This value 
of qc compares favorably with the estimate obtained via a real space renormalization group 
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method [|] (gc ~ 2.2), a method that yields the excellent value Qc = 4.08 in two dimensions; 
it is smaller than the estimates obtained from 1/q series expansions ^ [qc = 2.57±0.12 with 
a = 1.65±0.15 ) and from Monte-Carlo simulations g (g^ = 2.45±0.1 with u = 0.87±0.08). 
(Note however that the estimate of Ref. ^ is obtained by extrapolating numerical data in 
the range 2.8 < q < 3.2 .) In the case of the bcc lattice, the same value Qc ~ 2.14 (to within 
our numerical accuracy) is obtained from analysis of the latent heat and the correlation 
length, with exponents a ~ 1.61 and u ~ 1.11, respectively. On the other hand, no reliable 
results are obtained for the fee lattice. This is due probably to the fact that the asymptotic 
regime is not yet reached in the range 2.4 < q < 2.8, as can be seen in Fig. 4 with the small 
values of the correlation length ^. 

Further investigation of changeover behavior can be carried out by considering the q- 
dependence of the coordinates of the critical endpoint (restricting again the numerical study 
to g > 2.4). When q — » qc, one expects that he ^ and Tc{q) — Tec = Tc{qc) = Tt{qc)- 
As shown in Fig. 5, the data for the sc lattice are again well fitted by the simple form 
Tec - Tc{q) ~ hciqy with kBTcJiJc) = 0.775 ± 0.003 and p = 0.39 ± 0.02. (Note that the 
same power-law behavior with the mean- field exponent p = 0.4 = 2/5 is observed when a 
tricritical point is approached along a critical wing line [^.) As can be seen in Fig. 6, 
this result fits in with the global T — q phase diagram of the model on the sc lattice (this 
figure also includes the results for q < 2.4 for which no difference is found between Tt and 
Tc within the accuracy of our calculations). For T = Tec, we now find qc ~ 2.15, which is 
consistent with the preceding estimates. 

C. 1 < g < 2 and the bond-percolation limit 

When q < 2, the transition is clearly continuous, as illustrated in Fig. 7 for q = 2 and 
q = 1.8. One can thus investigate the asymptotic behavior of the zero-field susceptibility 
as the critical point is approached from the high-temperature phase. Fig. 8 shows log-log 
plots of the reduced susceptibility Xredil) = kBTx{q) / {q — 1) as a function of the reduced 
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temperature t = 1 — Tc{q)/T (the factor {q — 1) comes again from our definition of the 
magnetization). We also plot in the figure the results for q = 2.07 for which the susceptibility 
is apparently diverging, despite the fact that the transition is not really continuous, as 
discussed above and in Appendix B. In the whole range 1 < g < 2, the divergence of xIq) is 
governed asymptotically by the critical exponent 7 = 2 which is that of the spherical model. 
However, the results are affected by a strong crossover, and for t > 10^^ the divergence of 
x{q) is governed by an effective exponent that increases as q decreases (see Table III). For 
q = 2, this exponent is close to the exact Ising value 7 ~ 1.24, as noted in earlier work 
||T0| . Note that the increase of 7 as g decreases is the exact behavior of the two-dimensional 
model (see e.g, Ref. 0). 

An interesting limit is g ^ 1+, for which there is a well known mapping to the bond 



percolation problem ||2^: the critical temperature of the Potts model is related to the 
percolation treshold via the equation pc/ (1 —Pc) = e'^^"*-'^^^^ — 1 and 7 becomes the exponent 
for the mean cluster size. For q = 1.001, we find Pc = 0.241, 0.172 and 0.114 for the sc, bcc, 
and fee lattices, respectively. These values are in reasonable agreement with the current 
best estimates Pc = 0.2488,0.1803 and 0.119 [0. Moreover, in the range t > lO'^, the 
divergence of the susceptiblity is governed by the effective exponent 7 ~ 1.71, which is close 
to the expected value for three dimensional lattices, 7 ~ 1.80. 



V. CONCLUSION 

We have applied a thermodynamically self-consistent OZ approximation to the ferromag- 
netic Potts model on the various three-dimensional lattices.. The present theory gives an 
accurate description of the phase diagram of the model for all values of g > 1. The transition 
temperatures for instance are within 0.05% of the best known estimates for g = 3 and g = 4 
(first-order transitions) and even in the limit g ^ 1+ where the transition temperature gives 
the percolation threshold of the associated bond percolation problem (second-order transi- 
tion), they are within 4% of the Monte Carlo estimates for the sc, bcc, and fee lattices. The 
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SCOZA has an intrinsic limitation: the OZ assumption restricting the range of the direct 
correlation functions to that of the pair interaction potential prevents an accurate treatment 
of the asymptotic critical behavior. However, the region of concern appears to be very small 
(|T — Tc\/Tc < 10^^), which allows the theory to make good predictions for the effective 
critical exponents. One finds for instance 7 1.28 instead of 1.24 for q = 2 (Ising model) 
and 7 ~ 1.71 instead of 1.80 for g ^ 1+ (percolation). 

Another shortcoming that we have encountered in the present study results from the 
following fact: in order to describe the Potts model for arbitrary real values of q and thereby 
investigate the potential change in the order of the transition, one has to restrict the param- 
eter space to a subspace preserving the Zg_i symmetry between all but one of the q states. 
As a consequence of this restriction, the SCOZA approach becomes incomplete, with more 
unknowns than exact conditions to be satisfied. Additional approximations must thus be 
devised that go beyond the mere application of the OZ formalism. The approximation we 
chose in this work (see Eq. (43)) still allows, as stressed above, to locate with high accu- 
racy the transition temperatures and to provide a good description of the thermodynamic 
quantities over a wide range of the parameters (T, m, g). It also gives numerical indication 
that the order of the transition seems to change from second to first order at a qc slightly 
above 2. Yet, a proper description of this changeover phenomenon would require a better 
approximation than that used here. Making progress on this issue, that also plagues the 
SCOZA approach for systems with quenched disorder in the replica formalism (because of 
the necessity to make an analytic continuation for letting the number of replicas go to zero 
||T^), is certainly an important problem that must still be addressed to confirm the SCOZA 
as an accurate and a versatile theoretical scheme. 
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APPENDIX A: 



In this appendix, we derive the exact expressions of Cl^ and for m = q — 1. These 
expressions are used to build the approximate relationship between z and z*, Eq. (44). 

When m = g — 1, all spins are in state 1 and h = (3h oo. In order to get an expansion 
in powers of m = (g — 1 — m), we consider, at a given temperature, the excited states 
corresponding to configurations where one, two, or more spins are in a state different from 
1. The partition function is first expressed in powers oiy = e~^. To derive Eq. (^21) , we only 
need to consider the first and the second excited states. At this order, the partition function 
Z{f3J,h) in any dimension depends only on the coordination number c of the lattice. We 
find 



with X = e~^'^. We deduce from Eq. (|A1[) the expansion for the magnetization. 



= = 1)^"^ (A2) 

oh 

+q{q - l)x2'?^[(g - 1)(1 + c) + x-«(2 - q)c - ca;-^^]?/^? + o{y^'') . 
The inversion of Eq. ( |A2| ) gives the expansion of y'^ in powers of u, 

^ ~ qiq- + q'{q - l^x^^^ + ^^"^ ^^"^^ 

The expansion of the susceptibility x = dm/dh then follows, 

X = qu-^ ^ ^-u^ + Oiw^). (A4) 

q- 1 

Using keTx = G^\k = 0) = {q - 1 - m){l + m)[l + z + 0{z'^)] (see Eq. (39)), we then find 

. = 4^--' + fa- fa- 1)1, ^ o(u'). (A5) 
q{q-l) 

To derive Eq. (^2|), the expansion of z* is also needed. We then expand the enthalpy S, 
Eq. (^51), in powers of u, 
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where 2; = + O(u^), 2;* = + O(m^). Note that in writing Eq. (A6) as in deriving Eq. 
(A5), we have used the fact that the SCOZA structure with the direct correlation functions 
extending only to n.n. separation is exact at the order in u considered. Comparing Eq. (A6) 
to the expansion of £^ = —J d\n{Z)/d{pj)\f^ that is derived from Eqs. ( [Al| ) and (|A3|), 

we find 

^ x-'?[x-^-1]m + 0(m2) . (A8) 



q{q-l 
Since, from Eqs. (30,31), 

1 



and 



C^= [z^ + iq-2)zl] + 0iu) (A9) 

cq 



_2(5_J) , (AlO) 

cq 



we finally obtain 



Cl\\, m = g - 1) - Cf (A, m = g - 1) = ^(1 - e^'?/^)l (All) 
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APPENDIX B: 



In this appendix, we show that the proposed theory cannot truly predict a second-order 
transition (in zero external field) for q> 2 because of the approximation given in Eq. (43), 
that represents a linear interpolation formula for — Cp. 

Let us assume that the transition is continuous. Then, the inverse susceptibility = 
dh/dm for m = should go to zero at the transition. This implies that z 1 when m — 
at the critical temperature, and we may expand all quantities in the vicinity of z = 1 and 
m = 0. Since the 3-dimensional lattice Green's function behaves as 

P{z)r^Pm-ae + Oie% (Bl) 

where e = \/l — z and a is a positive constant [Q, we get from Eqs. (25), (39) and (40) 



that the enthalpy behaves as 
P(l)-1 



^ - ^T^ir^^ " (g-l)P(l) ^^ + " + ^^^^ 



where e* = \/\ — z* . 

As discussed in section IIIB, the zero- field condition at m = in the disordered phase 
implies that d£ / dm\\^rn=o = 0. From the above expression, this implies that 

d{e+{q~2)e*) 



dm 



= . (B3) 

A,m=0 



We must thus have 



e+(g-2)e* = 0(m^+") (B4) 

with u > 0. 

On the other hand, expanding Eq. (44) near the critical point yields 

aP{l) 

where AC(Ac), given by Eq. (^2|), is positive. 

There are two possibilities for Eqs. (B4) and (B5) to be compatible: 
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(1) Neither e nor e* have a term hnear in m, thereby implying via Eq.(B5) that P{1) — 
cAC(Ac) — 0. The numerical study, however, shows that for 2 < g < 2.4 this quantity is 
small but strictly positive in the range of temperatures where the transition takes place. 

(2) Both e and e* have terms linear in m, with Eq. (B5) satisfied (for q^2), and there 
is an exact cancellation of these terms so that Eq. (B4) is also obeyed. However, since e 
and e* are both non- negative quantities, this cancellation cannot occur for g > 2. 

Therefore, there cannot be a zero- field continuous transition for > 2. 
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TABLES 

TABLE I. Theoretical predictions for the inverse transition temperature, the dimensionless 

latent heat, the jump discontinuity in the order parameter, and the coordinates of the critical 
endpoint on the sc, bcc and fee lattices for the 3 and 4-state Potts models. 





lattice 


qPtJ 


Ae 


Am/{q - 1) 


qPcJ 




q=3 


sc 


0.5507 


0.199 


0.444 


0.5474 


0.00080 




bcc 


0.3957 


0.319 


0.452 


0.3926 


0.00135 




fee 


0.2576 


0.519 


0.457 


0.2556 


0.00124 


q=4 


sc 


0.6288 


0.576 


0.628 


0.6134 


0.00491 




bcc 


0.4549 


0.870 


0.636 


0.4418 


0.00622 




fee 


0.2970 


1.371 


0.638 


0.2880 


0.00670 
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TABLE II. Theoretical predictions for the inverse transition temperature, the dimensionless 
latent heat, the jump discontinuity in the order parameter, and the coordinates of the critical 
endpoint on the simple cubic lattice for 2.4 < q < 2.8 



Q 


qPtJ 


Ae 


Am/{q - 1) 


qPcJ 


Pchc 


2.8 


0.5321 


0.132 


0.387 


0.5302 


0.00039 


2.7 


0.5223 


0.099 


0.349 


0.5209 


0.00025 


2.6 


0.5122 


0.072 


0.314 


0.5113 


0.00014 


2.5 


0.5017 


0.045 


0.265 


0.5011 


0.00007 


2.4 


0.4908 


0.027 


0.226 


0.4904 


0.00003 
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TABLE III. Theoretical predictions for the effective critical exponent 7e//(9) that governs the 
divergence of the susceptibility for (T — Tc)/T > 10~^. 

9 leff_ 

2 1.28 

1.9 1.34 

1.8 1.39 

1.7 1.45 

1.001 1.71 
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FIGURES 




FIG. 1. Schematic plot of the spinodal curve (solid line) and the spontaneous magnetization 
versus temperature (dashed line) for second-order (a) and first-order (b) transitions. 
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FIG. 2. Spinodal curve (solid line) and spontaneous magnetization versus temperature (dashed 
line) for g = 3 (a) and q = 4 (b). 
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FIG. 4. Second-moment correlation length on the sc, bcc, and fee lattices for 2.4 < g < 3. 
The solid line is the algebraic fit, ^ ~ 4.61 {q — 2.14)^-^'^, of the data for the sc lattice 
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FIG. 5. Coordinates of the critical endpoint in the T — h plane for q = 2.4, 2.5, 2.6, 2.7, 2.8, 3, 
and 3.2 (from left to right). The corresponding log-log plot indicates that he vanishes at 
kBTcc/{Jc) = 0.775 ± 0.003. 
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FIG. 6. Predicted phase diagram of the g-state Potts model on the sc lattice. The solid 
and dashed curves represent the second-order and first-order lines, respectively. The second-order 
transition occurs in zero field only for q < Qc- (Note that within numerical accuracy, our results 
cannot distinguish first from second-order for 2 < g < 2.4.) 
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FIG. 7. Spinodal curve (solid line) and spontaneous magnetization versus temperature (dashed 
line) for g = 2 (a) and q = 1.8 (b). 
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iog,o(i-To/T) 



FIG. 8. Log-log plot of the reduced zero-field susceptibility Xred = ksTx/iq — 1) versus the 
reduced temperature t = \ — T^T for q = 2.07, 2, 1.9, 1.8, 1.7 (solid lines from bottom to top) and 
q = 1.001 (dashed line). 



36 



